%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: January 2021
%
% Purpose: Endogenous Targeting (step 1)
%   - Solve Initial equilibrium on nuc and nuv grids
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 

    % load estimates (exogenous targeting)
    load('../../output/estimates/fullmodel_estimates.mat', 'est');            
    
    
    
%% Initialize

    % N
    Nc = 19; %number of nu_c
    Nv = 11; %number of nu_v for each nu_c
       
    % nu_c
    nuc_matrix = kron(ones(1,Nv), linspace(1,0.1,Nc)');        
    
    % nu_v
    nuv_center = [3.050, 3.050, 3.050, 3.050, 3.050, 3.035, 3.025, 3.015, 3.000, 3.000,...
                  3.000, 2.985, 2.975, 2.950, 2.925, 2.900, 2.850, 2.780, 2.725];
    nuv_matrix = zeros(Nc,Nv);
    for i = 1:Nc
       nuv_matrix(i,:) = nuv_center(i) + 0.005*[-5:5];  
    end
         
    % placeholder
    fval_matrix = zeros(Nc,Nv);    
    converg.iter_inner = zeros(Nc,Nv);
    converg.iter_outer = zeros(Nc,Nv);
    converg.err_inner = zeros(Nc,Nv);
    converg.err_outer = zeros(Nc,Nv);
    converg.fix_id = cell(Nc,Nv);
    converg.fix_qindex = cell(Nc,Nv);
    converg.fix_qindex_history = cell(Nc,Nv);
    
    
%% Loop (Calibration)
for nc = 1:Nc    
for nv = 1:Nv
tic    

    nu_c = nuc_matrix(nc,nv);
    nu_v = nuv_matrix(nc,nv);
    fprintf('\nCalibration(%.0f,%.0f): nu_c = %.4f, nu_v = %.4f\n', nc, nv, nu_c, nu_v)
    
    %% Solve Initial Equilibrium

        % parameters
        define_param
        define_param_est_endogtarget     
        param.nu_c = nu_c; %calibrate nu_c
        param.nu_v = nu_v; %calibrate nu_v

        % data moment
        define_data_moment

        % constants
        define_omega_dist 
        define_const_endogtarget 

        % solve initial equilibrium
        eqm = solve_eqm_outer_endogtarget(param, const);    
        define_wage_schedule
        const.phi_v = eqm.phi_v; %add phi_v to const

        % calculate model moments and objfcn
        model = cal_model_moment(param, const, eqm); 
        bartik = cal_chg_bartik(param, const, eqm);     
        fval = cal_objfcn(data, model, bartik);



    %% Store results

        % display
        fprintf('param.nu_c = %.2f, param.nu_v = %.4f, fval = %.4f, n(qmax) = %.6f, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
            param.nu_c, param.nu_v, fval, eqm.nq(param.Q), eqm.err_outer, eqm.err_inner, eqm.iter_outer, eqm.iter_inner);    

        % store
        result.param = param;
        result.const = const; 
        result.eqm = eqm; 
        result.model = model; 
        result.bartik = bartik; 
        result.fval = fval; 

        % store
        fval_matrix(nc,nv)                  = fval;
        converg.iter_inner(nc,nv)           = eqm.iter_inner;
        converg.iter_outer(nc,nv)           = eqm.iter_outer;
        converg.err_inner(nc,nv)            = eqm.err_inner;
        converg.err_outer(nc,nv)            = eqm.err_outer;
        converg.fix_id{nc,nv}               = eqm.fix_id;
        converg.fix_qindex{nc,nv}           = eqm.fix_qindex;
        converg.fix_qindex_history{nc,nv}   = eqm.fix_qindex_history;        
        
toc        
end
end



%% Save

    mkdir('../../output/temp')    
    save('../../output/temp/endogtarget_initeqm_grid.mat', 'Nc','Nv','nuc_matrix', 'nuv_matrix', 'fval_matrix', 'converg')     


      